Analyse
Unpublished work
Baptiste Crinière-Boizet
library(tidyverse)
library(openxlsx)
library(furniture)
library(ggsci)
library(lmerTest)
library(emmeans)
library(interactions)
library(kableExtra)
library(forcats)
library(ggpubr)
library(FactoMineR)
library(factoextra)Context
This experiment was conducted on healthy subjects divided into two groups to assess the impact of a stressful situation on their food choices. The goal is to determine whether subjects exposed to stress tend to prefer junk food, an unhealthy but pleasurable diet.
The 146 subjects were divided into two groups: 74 in the stress group and 72 in the control group. The experiment took place at five time points, during which participants had to assess their stress levels and make food choices. At points 2 and 4, subjects in the stress group were subjected to a stressful experience.
Stress was assessed using four criteria: stress, shame, threat, and self-securement.
The first part of the study aims to assess whether the experiment produces the expected effects on participants’ stress.
Modelisation :
We model our problem using a mixed-effects model, including a random effect associated with the subject, as we have multiple measurements per subject, one at each time point.
\[ \text{Stress}_{ij} = \beta_0 + \alpha_{0i} + \beta_{1j} \cdot \text{time}_{ij} + \beta_{2i} \cdot \text{condition}_i + \beta_{3ij} \cdot \text{time}_{ij} \cdot \text{condition}_i + \epsilon_{ij} \]
Explanation of coefficients:
- \(\beta_0\): Constant term.
- \(\alpha_{0i}\): Individual random effect for subject \(i\)
- \(\beta_{1j}\): Fixed effect of the time variable
- \(\beta_{2i}\): Fixed effect of the condition for subject \(i\)
- \(\beta_{3ij}\): Interaction between the time variable and the condition
- \(\epsilon_{ij}\): Error term
# Importation des données
data <- openxlsx::read.xlsx("/home/baptiste.criniere/Documents/PB_FOOD_AH/Data/Cleaned_Data_SM2.xlsx", na.strings = c("NA"))
# Data management
data1 <- data %>%
dplyr::select(ID, FL_11_DO, VAS1_1, VAS2_1, VAS3_1, VAS4_1, VAS5_1) %>%
tidyr::pivot_longer(cols = c("VAS1_1", "VAS2_1", "VAS3_1", "VAS4_1", "VAS5_1"), names_to = "Time", values_to = "Stress") %>%
dplyr::mutate(Time = Time %>%
factor %>%
forcats::fct_recode("1" = "VAS1_1", "2" = "VAS2_1", "3" = "VAS3_1",
"4" = "VAS4_1", "5" = "VAS5_1"))Stress Study
data1 <- data %>%
dplyr::select(ID, FL_11_DO, VAS1_1, VAS2_1, VAS3_1, VAS4_1, VAS5_1) %>%
tidyr::pivot_longer(cols = c("VAS1_1", "VAS2_1", "VAS3_1", "VAS4_1", "VAS5_1"), names_to = "Time", values_to = "Stress") %>%
dplyr::mutate(Time = Time %>%
factor %>%
forcats::fct_recode("1" = "VAS1_1", "2" = "VAS2_1", "3" = "VAS3_1",
"4" = "VAS4_1", "5" = "VAS5_1"))Description
In this initial section, we describe the stress variable at each time point. In the following section, we will proceed with statistical inference to determine if these differences are solely due to ‘random chance.’
Table 1
In this first table, I have presented the average stress level for
each group at each time point, accompanied by the standard deviation in
parentheses.
We observe that the average stress level is higher for the stressed
group at each time point, including time 1. Additionally, we also note
that the variability in stress levels, represented by the standard
deviation, is higher for the stressed group.
tab <- data %>%
furniture::table1("Stress Time 1" = VAS1_1,
"Stress Time 2" = VAS2_1,
"Stress Time 3" = VAS3_1,
"Stress Time 4" = VAS4_1,
"Stress Time 5" = VAS5_1,
splitby =~ FL_11_DO,
na.rm = FALSE)knitr::kable(tab, caption = "") %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper()| . | ControlCondition | StressCondition |
|---|---|---|
| n = 72 | n = 74 | |
| Stress Time 1 | ||
| 13.9 (18.7) | 20.9 (27.1) | |
| Stress Time 2 | ||
| 12.0 (18.5) | 41.5 (35.4) | |
| Stress Time 3 | ||
| 12.2 (23.0) | 20.1 (26.4) | |
| Stress Time 4 | ||
| 24.0 (29.1) | 29.6 (34.9) | |
| Stress Time 5 | ||
| 18.3 (22.7) | 22.3 (30.6) |
Figure 1
In this figure, the density of the stress variable is represented for each condition and at each time.
data1 %>%
ggplot(aes(x = Stress))+
geom_density(aes(fill = FL_11_DO), alpha = 0.5)+
facet_grid(FL_11_DO ~ Time)+
theme_classic()+
scale_fill_jama()+
theme(legend.position = "none")
Figure 2
In this figure, we represent the stress variable at each time point using a boxplot. The two conditions (stress and control) are separated. The points correspond to the stress level of each individual and are connected over time by a line.
data1 %>%
ggplot(aes(x = Time, y = Stress))+
geom_point(size = 0.3, alpha = 0.25)+
geom_boxplot(aes(fill = Time), outlier.shape = NA, alpha = 0.65)+
geom_line(aes(group = ID), alpha = 0.25, size = 0.25)+
facet_wrap(~FL_11_DO)+
theme_bw()+
scale_fill_jama()+
theme(legend.position = "none")
Figure 3
This last figure is the same as the previous one but without the points and lines.
data1 %>%
ggplot(aes(x = Time, y = Stress))+
geom_boxplot(aes(fill = FL_11_DO), alpha = 0.5, outlier.shape = NA)+
theme_bw()+
scale_fill_jama()+
labs(fill = "Condition")
Modeling
In this section, we aim to determine whether the individual
differences observed in descriptive statistics (medians, means, figures)
are generalizable. To do this, we will use a mixed-effects model to test
the effect of time and condition (their interaction) on the stress
score.
The equation of the model is as follows: Stress = Time*Condition +
(1|ID)
Global Test
First, we test the overall effect of each variable. The p-values associated with the overall significance of each coefficient are below 0.05, the significance threshold. This indicates that there is an effect of time, condition, and their interaction on the stress score. You can read the p-value associated with the significance test of each variable in the last column of the following table.
model <- lmerTest::lmer(Stress ~ Time*FL_11_DO + (1|ID), data = data1)
# summary(model)
anova(model)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Time 14632.0 3658.0 4 566.22 12.6562 7.168e-10 ***
## FL_11_DO 2408.6 2408.6 1 143.72 8.3335 0.004493 **
## Time:FL_11_DO 16254.5 4063.6 4 566.22 14.0596 6.048e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Comparison of the Two Groups at Each Time Point
In this section, our goal is to compare the two groups at each time point. We observe only one significant difference between the two groups, which occurs at time T2. At that moment, the “stressed” group has an estimated mean value of the stress score higher by 29.57 compared to the “control” group (p-value < 0.0001).
emm <- emmeans(model, ~ FL_11_DO | Time)
contrast(emm, interaction = c("pairwise"))## Time = 1:
## FL_11_DO_pairwise estimate SE df t.ratio p.value
## ControlCondition - StressCondition -6.93 4.53 286 -1.530 0.1272
##
## Time = 2:
## FL_11_DO_pairwise estimate SE df t.ratio p.value
## ControlCondition - StressCondition -29.57 4.53 286 -6.523 <.0001
##
## Time = 3:
## FL_11_DO_pairwise estimate SE df t.ratio p.value
## ControlCondition - StressCondition -8.31 4.61 302 -1.803 0.0724
##
## Time = 4:
## FL_11_DO_pairwise estimate SE df t.ratio p.value
## ControlCondition - StressCondition -5.62 4.53 286 -1.240 0.2159
##
## Time = 5:
## FL_11_DO_pairwise estimate SE df t.ratio p.value
## ControlCondition - StressCondition -4.03 4.53 286 -0.890 0.3744
##
## Degrees-of-freedom method: kenward-rogerComparison of Times for Each Group
In this section, the goal is to compare the stress score values between different time points for each group. I conducted comparisons between stress levels at one time compared to the previous time.
We observe that there is a significant difference for the control group at time 4 compared to times 3 and 5.
Furthermore, we note that there is a significant difference at each time point for the “stressed” group.
emm2 <- emmeans(model, ~ Time | FL_11_DO)
contrasts <- list("T1 vs T2" = c(1, -1, 0, 0, 0),
"T2 vs T3" = c(0, 1, -1, 0, 0),
"T3 vs T4" = c(0, 0, 1, -1, 0),
"T4 vs T5" = c(0, 0, 0, 1, -1))
contrast(emm2, contrasts)## FL_11_DO = ControlCondition:
## contrast estimate SE df t.ratio p.value
## T1 vs T2 1.98611 2.83 566 0.701 0.4836
## T2 vs T3 -0.00571 2.88 567 -0.002 0.9984
## T3 vs T4 -12.02207 2.88 567 -4.169 <.0001
## T4 vs T5 5.69444 2.83 566 2.010 0.0449
##
## FL_11_DO = StressCondition:
## contrast estimate SE df t.ratio p.value
## T1 vs T2 -20.64865 2.79 566 -7.388 <.0001
## T2 vs T3 21.25137 2.87 568 7.407 <.0001
## T3 vs T4 -9.33245 2.87 568 -3.253 0.0012
## T4 vs T5 7.28378 2.79 566 2.606 0.0094
##
## Degrees-of-freedom method: kenward-rogerFinal Figure
Here is the final figure representing the estimated marginal value of the stress score by condition and time with a confidence interval.
interactions::cat_plot(model,
pred = "Time",
modx = "FL_11_DO",
interval = TRUE,
y.label = "Estimated Marginal Mean Stress Score",
x.label = "Time",
legend.main = "Condition:",
modx.labels = c("Control", "Stress"),
colors = c("gray40", "black"),
point.shape = TRUE,
pred.point.size = 5,
dodge.width = .3,
errorbar.width = .25,
geom = "line") +
theme_bw() +
theme(legend.key.width = unit(2, "cm"),
legend.background = element_rect(color = "Black"),
legend.position = c(1, 0.72),
legend.justification = c(1.1, -0.1)) +
scale_y_continuous(breaks = seq(from = 0, to = 50, by = 10))+
labs(title = "Modelisation results")
Shame Study
In this first part, we describe the variable “ashamed” at each time point. In the following section, we will proceed with statistical inference to determine if these differences are solely due to “random chance.”
Description
data1 <- data %>%
dplyr::select(ID, FL_11_DO, VAS1_4, VAS2_4, VAS3_4, VAS4_4, VAS5_4) %>%
tidyr::pivot_longer(cols = c("VAS1_4", "VAS2_4", "VAS3_4", "VAS4_4", "VAS5_4"), names_to = "Time", values_to = "Ashamed") %>%
dplyr::mutate(Time = Time %>%
factor %>%
forcats::fct_recode("1" = "VAS1_4", "2" = "VAS2_4", "3" = "VAS3_4",
"4" = "VAS4_4", "5" = "VAS5_4"))Table 1
In this first table, I have represented the average “ashamed” value for each group at each time point, accompanied by the standard deviation in parentheses. We observe that the average “ashamed” value is higher for the stressed group at times 1, 2, 3, 4, except at time 5.
tab <- data %>%
furniture::table1("Ashamed Time 1" = VAS1_4,
"Ashamed Time 2" = VAS2_4,
"Ashamed Time 3" = VAS3_4,
"Ashamed Time 4" = VAS4_4,
"Ashamed Time 5" = VAS5_4,
splitby =~ FL_11_DO,
na.rm = FALSE)knitr::kable(tab, caption = "") %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper()| . | ControlCondition | StressCondition |
|---|---|---|
| n = 72 | n = 74 | |
| Ashamed Time 1 | ||
| 4.8 (16.4) | 6.3 (18.3) | |
| Ashamed Time 2 | ||
| 5.2 (16.9) | 17.3 (24.9) | |
| Ashamed Time 3 | ||
| 4.8 (17.1) | 6.1 (16.8) | |
| Ashamed Time 4 | ||
| 11.1 (24.6) | 13.3 (25.9) | |
| Ashamed Time 5 | ||
| 8.6 (19.5) | 8.2 (19.4) |
Figure 1
In this figure, we represent the “ashamed” variable at each time point using a boxplot. The two conditions (stress and control) are separated. The points correspond to the level of “ashamed” for each individual and are connected over time by a line.
data1 %>%
ggplot(aes(x = Time, y = Ashamed))+
geom_point(size = 0.3, alpha = 0.25)+
geom_boxplot(aes(fill = Time), outlier.shape = NA, alpha = 0.65)+
geom_line(aes(group = ID), alpha = 0.25, size = 0.25)+
facet_wrap(~FL_11_DO)+
theme_bw()+
scale_fill_jama()+
theme(legend.position = "none")
Figure 2
data1 %>%
ggplot(aes(x = Time, y = Ashamed))+
geom_boxplot(aes(fill = FL_11_DO), alpha = 0.5, outlier.shape = NA)+
theme_bw()+
scale_fill_jama()+
labs(fill = "Condition")
Modeling
The equation of the model is as follows: Ashamed = Time*Condition + (1|ID)
Global Test
Firstly, we test the overall effect of each variable. The p-values associated with the global significance of time and the time-condition interaction are below 0.05, the significance threshold. This indicates that there is an effect of time and the interaction on the “ashamed” variable.
model <- lmerTest::lmer(Ashamed ~ Time*FL_11_DO + (1|ID), data = data1)
# summary(model)
anova(model)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Time 5187.9 1296.98 4 566.07 9.1021 3.941e-07 ***
## FL_11_DO 229.8 229.79 1 143.65 1.6127 0.2062
## Time:FL_11_DO 3489.5 872.36 4 566.07 6.1222 7.928e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Comparison of the Two Groups at Each Time Point
In this section, our objective is to compare the two groups at each time point. We observe only one significant difference between the two groups, which occurs at time T2 (p-value < 0.0001).
emm <- emmeans(model, ~ FL_11_DO | Time)
contrast(emm, interaction = c("pairwise"))## Time = 1:
## FL_11_DO_pairwise estimate SE df t.ratio p.value
## ControlCondition - StressCondition -1.574 3.38 262 -0.466 0.6419
##
## Time = 2:
## FL_11_DO_pairwise estimate SE df t.ratio p.value
## ControlCondition - StressCondition -12.145 3.38 262 -3.592 0.0004
##
## Time = 3:
## FL_11_DO_pairwise estimate SE df t.ratio p.value
## ControlCondition - StressCondition -2.807 3.43 276 -0.818 0.4142
##
## Time = 4:
## FL_11_DO_pairwise estimate SE df t.ratio p.value
## ControlCondition - StressCondition -2.187 3.38 262 -0.647 0.5182
##
## Time = 5:
## FL_11_DO_pairwise estimate SE df t.ratio p.value
## ControlCondition - StressCondition 0.395 3.38 262 0.117 0.9070
##
## Degrees-of-freedom method: kenward-rogerComparison of Time Points for Each Group
We observe a significant difference for the control group at time 4 compared to times 3. Additionally, we notice a significant difference at each time point for the “stressed” group.
emm2 <- emmeans(model, ~ Time | FL_11_DO)
contrasts <- list("T1 vs T2" = c(1, -1, 0, 0, 0),
"T2 vs T3" = c(0, 1, -1, 0, 0),
"T3 vs T4" = c(0, 0, 1, -1, 0),
"T4 vs T5" = c(0, 0, 0, 1, -1))
contrast(emm2, contrasts)## FL_11_DO = ControlCondition:
## contrast estimate SE df t.ratio p.value
## T1 vs T2 -0.389 1.99 566 -0.195 0.8451
## T2 vs T3 0.653 2.02 567 0.323 0.7472
## T3 vs T4 -6.570 2.02 567 -3.244 0.0012
## T4 vs T5 2.444 1.99 566 1.229 0.2197
##
## FL_11_DO = StressCondition:
## contrast estimate SE df t.ratio p.value
## T1 vs T2 -10.959 1.96 566 -5.585 <.0001
## T2 vs T3 9.991 2.01 567 4.959 <.0001
## T3 vs T4 -5.950 2.01 567 -2.954 0.0033
## T4 vs T5 5.027 1.96 566 2.562 0.0107
##
## Degrees-of-freedom method: kenward-rogerFinal Figure
Here is the final figure representing the estimated marginal value of “ashamed” by condition and by time with a confidence interval.
interactions::cat_plot(model,
pred = "Time",
modx = "FL_11_DO",
interval = TRUE,
y.label = "Estimated Marginal Mean Ashamed Score",
x.label = "Time",
legend.main = "Condition:",
modx.labels = c("Control", "Stress"),
colors = c("gray40", "black"),
point.shape = TRUE,
pred.point.size = 5,
dodge.width = .3,
errorbar.width = .25,
geom = "line") +
theme_bw() +
theme(legend.key.width = unit(2, "cm"),
legend.background = element_rect(color = "Black"),
legend.position = c(1, 0.72),
legend.justification = c(1.1, -0.1)) +
scale_y_continuous(breaks = seq(from = 0, to = 50, by = 10))+
labs(title = "Modelisation results")
Study of Threatened
We describe the Threatened variable at each time point.
Description
data1 <- data %>%
dplyr::select(ID, FL_11_DO, VAS1_5, VAS2_5, VAS3_5, VAS4_5, VAS5_5) %>%
tidyr::pivot_longer(cols = c("VAS1_5", "VAS2_5", "VAS3_5", "VAS4_5", "VAS5_5"), names_to = "Time", values_to = "Threatened") %>%
dplyr::mutate(Time = Time %>%
factor %>%
forcats::fct_recode("1" = "VAS1_5", "2" = "VAS2_5", "3" = "VAS3_5",
"4" = "VAS4_5", "5" = "VAS5_5"))Table 1
In this first table, I have presented the mean value of Threatened for each group at each time point, accompanied by the standard deviation in parentheses. It can be observed that the mean of Threatened is higher for the stressed group at each time point, including at time 1.
tab <- data %>%
furniture::table1("Threatened Time 1" = VAS1_5,
"Threatened Time 2" = VAS2_5,
"Threatened Time 3" = VAS3_5,
"Threatened Time 4" = VAS4_5,
"Threatened Time 5" = VAS5_5,
splitby =~ FL_11_DO,
na.rm = FALSE)knitr::kable(tab, caption = "") %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper()| . | ControlCondition | StressCondition |
|---|---|---|
| n = 72 | n = 74 | |
| Threatened Time 1 | ||
| 2.5 (10.0) | 4.3 (13.0) | |
| Threatened Time 2 | ||
| 2.8 (11.1) | 19.0 (27.9) | |
| Threatened Time 3 | ||
| 2.5 (10.3) | 5.8 (16.2) | |
| Threatened Time 4 | ||
| 3.2 (7.4) | 9.4 (21.2) | |
| Threatened Time 5 | ||
| 4.2 (9.6) | 6.1 (16.3) |
Figure 1
data1 %>%
ggplot(aes(x = Time, y = Threatened))+
geom_point(size = 0.3, alpha = 0.25)+
geom_boxplot(aes(fill = Time), outlier.shape = NA, alpha = 0.65)+
geom_line(aes(group = ID), alpha = 0.25, size = 0.25)+
facet_wrap(~FL_11_DO)+
theme_bw()+
scale_fill_jama()+
theme(legend.position = "none")
Figure 2
data1 %>%
ggplot(aes(x = Time, y = Threatened))+
geom_boxplot(aes(fill = FL_11_DO), alpha = 0.5, outlier.shape = NA)+
theme_bw()+
scale_fill_jama()+
labs(fill = "Condition")
Modeling
The equation of the model is as follows: Threatened = Time*Condition + (1|ID)
Global Test
Initially, we test the overall effect of each variable. The p-values associated with the overall significance of each coefficient are below 0.05, the significance threshold. This indicates that there is an effect of time, condition, and their interaction on Threatened.
model <- lmerTest::lmer(Threatened ~ Time*FL_11_DO + (1|ID), data = data1)
# summary(model)
anova(model)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Time 5008.5 1252.11 4 566.78 11.8571 2.943e-09 ***
## FL_11_DO 858.2 858.19 1 144.18 8.1268 0.005002 **
## Time:FL_11_DO 5278.7 1319.67 4 566.78 12.4969 9.489e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Comparison of Two Groups at Each Time Point
We observe two statistically significant differences between the two groups, occurring at T2 and T4.
emm <- emmeans(model, ~ FL_11_DO | Time)
contrast(emm, interaction = c("pairwise"))## Time = 1:
## FL_11_DO_pairwise estimate SE df t.ratio p.value
## ControlCondition - StressCondition -1.82 2.57 316 -0.709 0.4791
##
## Time = 2:
## FL_11_DO_pairwise estimate SE df t.ratio p.value
## ControlCondition - StressCondition -16.22 2.57 316 -6.300 <.0001
##
## Time = 3:
## FL_11_DO_pairwise estimate SE df t.ratio p.value
## ControlCondition - StressCondition -3.57 2.62 334 -1.359 0.1750
##
## Time = 4:
## FL_11_DO_pairwise estimate SE df t.ratio p.value
## ControlCondition - StressCondition -6.18 2.57 316 -2.402 0.0169
##
## Time = 5:
## FL_11_DO_pairwise estimate SE df t.ratio p.value
## ControlCondition - StressCondition -1.85 2.57 316 -0.717 0.4740
##
## Degrees-of-freedom method: kenward-rogerComparison of Times for Each Group
For the stressed group, there are only two significant differences: between time 1 and 2, and between time 2 and 3.
emm2 <- emmeans(model, ~ Time | FL_11_DO)
contrasts <- list("T1 vs T2" = c(1, -1, 0, 0, 0),
"T2 vs T3" = c(0, 1, -1, 0, 0),
"T3 vs T4" = c(0, 0, 1, -1, 0),
"T4 vs T5" = c(0, 0, 0, 1, -1))
contrast(emm2, contrasts)## FL_11_DO = ControlCondition:
## contrast estimate SE df t.ratio p.value
## T1 vs T2 -0.292 1.71 566 -0.170 0.8648
## T2 vs T3 0.308 1.74 567 0.177 0.8597
## T3 vs T4 -0.683 1.74 567 -0.392 0.6952
## T4 vs T5 -1.042 1.71 566 -0.608 0.5433
##
## FL_11_DO = StressCondition:
## contrast estimate SE df t.ratio p.value
## T1 vs T2 -14.689 1.69 566 -8.695 <.0001
## T2 vs T3 12.963 1.73 568 7.476 <.0001
## T3 vs T4 -3.301 1.73 568 -1.904 0.0574
## T4 vs T5 3.297 1.69 566 1.952 0.0515
##
## Degrees-of-freedom method: kenward-rogerFinal Figure
Here is the final figure representing the estimated marginal value of Threatened by condition and time with a confidence interval.
interactions::cat_plot(model,
pred = "Time",
modx = "FL_11_DO",
interval = TRUE,
y.label = "Estimated Marginal Mean Threatened Score",
x.label = "Time",
legend.main = "Condition:",
modx.labels = c("Control", "Stress"),
colors = c("gray40", "black"),
point.shape = TRUE,
pred.point.size = 5,
dodge.width = .3,
errorbar.width = .25,
geom = "line") +
theme_bw() +
theme(legend.key.width = unit(2, "cm"),
legend.background = element_rect(color = "Black"),
legend.position = c(1, 0.72),
legend.justification = c(1.1, -0.1)) +
scale_y_continuous(breaks = seq(from = 0, to = 50, by = 10))+
labs(title = "Modelisation results")
Study of Self-secured
We describe the variable Self-secured at each point in time.
Description
data1 <- data %>%
dplyr::select(ID, FL_11_DO, VAS1_6, VAS2_6, VAS3_6, VAS4_6, VAS5_6) %>%
tidyr::pivot_longer(cols = c("VAS1_6", "VAS2_6", "VAS3_6", "VAS4_6", "VAS5_6"), names_to = "Time", values_to = "Self_secured") %>%
dplyr::mutate(Time = Time %>%
factor %>%
forcats::fct_recode("1" = "VAS1_6", "2" = "VAS2_6", "3" = "VAS3_6",
"4" = "VAS4_6", "5" = "VAS5_6"))Table 1
In this first table, I have presented the average value of Self-secured for each group at each point in time, accompanied by the standard deviation in parentheses. It is observed that the average of Self-secured is lower for the stressed group at each point in time, including time 1.
tab <- data %>%
furniture::table1("Self secured Time 1" = VAS1_6,
"Self secured Time 2" = VAS2_6,
"Self secured Time 3" = VAS3_6,
"Self secured Time 4" = VAS4_6,
"Self secured Time 5" = VAS5_6,
splitby =~ FL_11_DO,
na.rm = FALSE)knitr::kable(tab, caption = "") %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper()| . | ControlCondition | StressCondition |
|---|---|---|
| n = 72 | n = 74 | |
| Self secured Time 1 | ||
| 68.4 (27.3) | 63.1 (33.0) | |
| Self secured Time 2 | ||
| 69.8 (27.2) | 46.7 (36.3) | |
| Self secured Time 3 | ||
| 69.0 (30.0) | 58.4 (33.8) | |
| Self secured Time 4 | ||
| 62.1 (31.2) | 53.0 (36.7) | |
| Self secured Time 5 | ||
| 60.9 (30.0) | 57.0 (34.3) |
Figure 1
data1 %>%
ggplot(aes(x = Time, y = Self_secured))+
geom_point(size = 0.3, alpha = 0.25)+
geom_boxplot(aes(fill = Time), outlier.shape = NA, alpha = 0.65)+
geom_line(aes(group = ID), alpha = 0.25, size = 0.25)+
facet_wrap(~FL_11_DO)+
theme_bw()+
scale_fill_jama()+
theme(legend.position = "none")
Figure 2
data1 %>%
ggplot(aes(x = Time, y = Self_secured))+
geom_boxplot(aes(fill = FL_11_DO), alpha = 0.5, outlier.shape = NA)+
theme_bw()+
scale_fill_jama()+
labs(fill = "Condition")
Modeling
The equation of the model is as follows: Self-secured = Time*Condition + (1|ID)
Global Test
Initially, we test the overall effect of each variable. The p-values associated with the overall significance of each coefficient are less than 0.05, the significance threshold. This indicates that there is an effect of time, condition, and their interaction on Self-secured.
model <- lmerTest::lmer(Self_secured ~ Time*FL_11_DO + (1|ID), data = data1)
# summary(model)
anova(model)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Time 7960.5 1990.13 4 566.15 8.5222 1.108e-06 ***
## FL_11_DO 998.4 998.38 1 143.91 4.2753 0.04046 *
## Time:FL_11_DO 8450.1 2112.54 4 566.15 9.0464 4.352e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Comparison of the Two Groups at Each Time Point
We observe two statistically significant differences between the two groups that occur at T2.
emm <- emmeans(model, ~ FL_11_DO | Time)
contrast(emm, interaction = c("pairwise"))## Time = 1:
## FL_11_DO_pairwise estimate SE df t.ratio p.value
## ControlCondition - StressCondition 5.36 5.34 211 1.005 0.3162
##
## Time = 2:
## FL_11_DO_pairwise estimate SE df t.ratio p.value
## ControlCondition - StressCondition 23.08 5.34 211 4.323 <.0001
##
## Time = 3:
## FL_11_DO_pairwise estimate SE df t.ratio p.value
## ControlCondition - StressCondition 8.52 5.39 219 1.579 0.1157
##
## Time = 4:
## FL_11_DO_pairwise estimate SE df t.ratio p.value
## ControlCondition - StressCondition 9.07 5.34 211 1.699 0.0908
##
## Time = 5:
## FL_11_DO_pairwise estimate SE df t.ratio p.value
## ControlCondition - StressCondition 3.98 5.34 211 0.747 0.4562
##
## Degrees-of-freedom method: kenward-rogerComparison of Times for Each Group
For the stressed group, there are two significant differences at almost all times except for the comparison between times 4 and 5. Meanwhile, for the control group, there is a difference between times 3 and 4.
emm2 <- emmeans(model, ~ Time | FL_11_DO)
contrasts <- list("T1 vs T2" = c(1, -1, 0, 0, 0),
"T2 vs T3" = c(0, 1, -1, 0, 0),
"T3 vs T4" = c(0, 0, 1, -1, 0),
"T4 vs T5" = c(0, 0, 0, 1, -1))
contrast(emm2, contrasts)## FL_11_DO = ControlCondition:
## contrast estimate SE df t.ratio p.value
## T1 vs T2 -1.40 2.55 566 -0.551 0.5820
## T2 vs T3 1.51 2.59 567 0.583 0.5600
## T3 vs T4 6.25 2.59 567 2.411 0.0162
## T4 vs T5 1.11 2.55 566 0.436 0.6628
##
## FL_11_DO = StressCondition:
## contrast estimate SE df t.ratio p.value
## T1 vs T2 16.31 2.51 566 6.493 <.0001
## T2 vs T3 -13.05 2.58 567 -5.058 <.0001
## T3 vs T4 6.80 2.58 567 2.638 0.0086
## T4 vs T5 -3.97 2.51 566 -1.581 0.1143
##
## Degrees-of-freedom method: kenward-rogerFinal Figure
Here is the final figure representing the estimated marginal value of Self-secured by condition and over time with a confidence interval.
interactions::cat_plot(model,
pred = "Time",
modx = "FL_11_DO",
interval = TRUE,
y.label = "Estimated Marginal Mean Self-secure Score",
x.label = "Time",
legend.main = "Condition:",
modx.labels = c("Control", "Stress"),
colors = c("gray40", "black"),
point.shape = TRUE,
pred.point.size = 5,
dodge.width = .3,
errorbar.width = .25,
geom = "line") +
theme_bw() +
theme(legend.key.width = unit(2, "cm"),
legend.background = element_rect(color = "Black"),
legend.position = c(1, 0.72),
legend.justification = c(1.1, -0.1)) +
scale_y_continuous(breaks = seq(from = 0, to = 50, by = 10))+
labs(title = "Modelisation results")
Correlation
I examined the correlation at each time point between the different variables.
T1
data1 <- data %>%
dplyr::select(VAS1_1, VAS1_4, VAS1_5, VAS1_6) %>%
dplyr::rename(Stress = VAS1_1,
Ashamed = VAS1_4,
Threatened = VAS1_5,
'Self secured' = VAS1_6,
)
M <- cor(data1)
corrplot::corrplot(M, type = "upper", order = "original",
tl.col = "black", tl.srt = 45, addCoef.col = "White")
T2
data1 <- data %>%
dplyr::select(VAS2_1, VAS2_4, VAS2_5, VAS2_6) %>%
dplyr::rename(Stress = VAS2_1,
Ashamed = VAS2_4,
Threatened = VAS2_5,
'Self secured' = VAS2_6,
)
M <- cor(data1)
corrplot::corrplot(M, type = "upper", order = "original",
tl.col = "black", tl.srt = 45, addCoef.col = "White")
T3
data1 <- data %>%
dplyr::select(VAS3_1, VAS3_4, VAS3_5, VAS3_6) %>%
dplyr::rename(Stress = VAS3_1,
Ashamed = VAS3_4,
Threatened = VAS3_5,
'Self secured' = VAS3_6,
) %>%
dplyr::filter(if_any(everything(), ~!is.na(.)))
M <- cor(data1)
corrplot::corrplot(M, type = "upper", order = "original",
tl.col = "black", tl.srt = 45, addCoef.col = "White")
T4
data1 <- data %>%
dplyr::select(VAS4_1, VAS4_4, VAS4_5, VAS4_6) %>%
dplyr::rename(Stress = VAS4_1,
Ashamed = VAS4_4,
Threatened = VAS4_5,
'Self secured' = VAS4_6,
) %>%
dplyr::filter(if_any(everything(), ~!is.na(.)))
M <- cor(data1)
corrplot::corrplot(M, type = "upper", order = "original",
tl.col = "black", tl.srt = 45, addCoef.col = "White")
T5
data1 <- data %>%
dplyr::select(VAS5_1, VAS5_4, VAS5_5, VAS5_6) %>%
dplyr::rename(Stress = VAS5_1,
Ashamed = VAS5_4,
Threatened = VAS5_5,
'Self secured' = VAS5_6,
) %>%
dplyr::filter(if_any(everything(), ~!is.na(.)))
M <- cor(data1)
corrplot::corrplot(M, type = "upper", order = "original",
tl.col = "black", tl.srt = 45, addCoef.col = "White")
PCA: Principal Component Analysis
We observe that overall, for each time point, the variables related to the questions “Hostile,” “Upset,” “Ashamed,” “Nervous,” “Afraid” group together, indicating that they are correlated. Meanwhile, the variables “Alert,” “Inspired,” “Determined,” “Active,” “Attentive” group together, indicating their correlation. Additionally, these two groups of variables are orthogonal, suggesting independence between these two sets of questions.
T1
mat <- data %>%
dplyr::select(starts_with("PANAS1"), FL_11_DO)
names(mat) <- c("Upset", "Hostile", "Alert", "Ashamed", "Inspired", "Nervous", "Determined", "Attentive", "Afraid", "Active", "Group")
res.pca <- PCA(mat[,-11], scale.unit = TRUE, graph = FALSE)
fig1 <- fviz_pca_var(res.pca, col.var = "black")
fig1
T2
mat <- data %>%
dplyr::select(starts_with("PANAS2"), FL_11_DO)
names(mat) <- c("Upset", "Hostile", "Alert", "Ashamed", "Inspired", "Nervous", "Determined", "Attentive", "Afraid", "Active", "Group")
res.pca <- PCA(mat[,-11], scale.unit = TRUE, graph = FALSE)
fig1 <- fviz_pca_var(res.pca, col.var = "black")
fig1
T3
mat <- data %>%
dplyr::select(starts_with("PANAS3"))
names(mat) <- c("Upset", "Hostile", "Alert", "Ashamed", "Inspired", "Nervous", "Determined", "Attentive", "Afraid", "Active")
res.pca <- PCA(mat, scale.unit = TRUE, graph = FALSE)
fig1 <- fviz_pca_var(res.pca, col.var = "black")
fig1
T4
mat <- data %>%
dplyr::select(starts_with("PANAS4"))
names(mat) <- c("Upset", "Hostile", "Alert", "Ashamed", "Inspired", "Nervous", "Determined", "Attentive", "Afraid", "Active")
res.pca <- PCA(mat, scale.unit = TRUE, graph = FALSE)
fig1 <- fviz_pca_var(res.pca, col.var = "black")
fig1
T5
mat <- data %>%
dplyr::select(starts_with("PANAS5"))
names(mat) <- c("Upset", "Hostile", "Alert", "Ashamed", "Inspired", "Nervous", "Determined", "Attentive", "Afraid", "Active")
res.pca <- PCA(mat, scale.unit = TRUE, graph = FALSE)
fig1 <- fviz_pca_var(res.pca, col.var = "black")
fig1
MFA : multiple factorial analysis
mat <- data %>%
dplyr::select(starts_with("PANAS")) %>%
dplyr::select(matches("_1$"), matches("_2$"), matches("_3$"), matches("_4$"), matches("_5$"),
matches("_6$"), matches("_7$"), matches("_8$"), matches("_9$"), matches("_10$"))
MFA(mat, group = c(5, 5, 5, 5, 5, 5, 5, 5, 5, 5), type = rep("c", 10),
name.group = c("Upset", "Hostile", "Alert", "Ashamed", "Inspired", "Nervous", "Determined", "Attentive", "Afraid", "Active"))
